function [ct, lags,na, nb, dt] = simxcorr2(varargin)
%
% simxcorr2 - calculate crosscorrelation of input arrays, averaging across
% multiple trials/observations. 
% input data is event time lists, as either simple array or as a cell
% array, e.g.:
% a=[1:100];
% b = 1+1;
% simxcorr2(a, b);
% or:
% a = {[1:100]; [2:3:150]};
% b = a+1;
% simxcorr2(a, b);
% Autocorrelation is handled as 
global CONTROL
tmaxlag = 2000;
if(nargin == 0)
    sf = getmainselection;

    A=CONTROL(sf).EPSC1;
    B=CONTROL(sf).EPSC3;
    D1 = A.eventtime;
    D2 = B.eventtime;
else
    D1 = varargin{1};
    D2 = varargin{2};
    if(iscell(D1)) % assume both d1 and d2 are cell arrays
        maxlat = max(max(cell2mat(D1)), max(cell2mat(D2)));
    else
        maxlat = max(max(D1), max(D2));
    end;
    
end;
tmaxlag = maxlat;
if(nargin >= 3 & ishandle(varargin{3}) & ~isempty(varargin{3})) % third argument would be handle to an axes for the plot
    haxes = varargin{3};
else
    haxes = [];
end;
minlat = 0;
if(nargin >= 4 & length(varargin{4} == 2))
    awin = varargin{4};
    minlat = min(awin);
    maxlat = max(awin);
end;
if(nargin >= 5)
    tmaxlag = varargin{5};
end;
%nargin
c={}; lags = {};
dt = 1;
twin = tmaxlag;

maxlag = twin; % floor(twin/dt);
tmax = 0;
atmax = 0;
btmax = 0;
for i = 1:size(D1, 2)
    if(iscell(D1))
        t1 = max([D1{i}]/dt);
        t2 = max([D2{i}]/dt);
    else
        t1 = max(D1/dt);
        t2 = max(D2/dt);
    end;
    if(t1 > atmax)
        atmax = t1;
    end;
    if(t2 > btmax)
        btmax = t2;
    end;
end;
if(atmax > btmax)
    tmax = atmax;
else
    tmax = btmax;
end;

zt = zeros(1, floor(tmax)+1); % build sample base.
ct = zeros(1, maxlag*2+1);
na = 0;
nb = 0;
list = [2:size(D1,2)];
if(isempty(list))
    list = [size(D1,2)];
end;
for i = list
    if(iscell(D1))
        a = floor([D1{i}]/dt);
    else
        a = floor(D1/dt);
    end;
    a = a(find(a > minlat & a < maxlat));
    if(iscell(D2))
        b = floor([D2{i}]/dt);
    else
        b = floor(D2/dt);
    end;
    b = b(find(b > minlat & b < maxlat));
    ia = zt;
    ia(a) = 1;
    ib = zt;
    ib(b) = 1;

    if(any(ia ~= 0) & any (ib ~= 0))
        na = na + length(find(ia == 1));
        nb = nb + length(find(ib == 1));
        [cb, lags] = xcorr(ia, ib, maxlag, 'none');
        ct = ct + cb;
    end;


end;
% ct = ct/size(D1,2); % average....
ct = 1000*ct/(dt*na); % normalizing by input event count ... to get rate out (in spikes per second).
if(isempty(haxes))
 %   return;
    h = findobj('tag', 'simcorrfig');
    if(isempty(h))
        h = figure;
        set(h, 'tag', 'simcorrfig');
    else
        figure(h);
        clf;
    end;
else
    axes(haxes); % select the axes handle
end

%stem(lags, ct, 'filled', 'linewidth', 2.0, 'markersize', 2.0, 'marker', 's');
bar(lags, ct);
set(gca, 'Xlim', [-tmaxlag tmaxlag]);
%set(gca, 'Ylim', [0 0.2]);


